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SUMMARY 

A  computer  programme  is  described  by  means  of  which  one  may  estimate 
the  potential  accuracy  of  the  elements  of  a  given  satellite  orbit,  if 
determined  from  observational  data  of  specified  type  and  assumed  accuracy. 
An  application  of  the  programme  is  made  to  an  orbit  of  six  hours  period 
determined  from  radar  observations  at  a  single  station. 
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1  INTRODUCTION 

There  is  a  recurrent  need,  in  assessment  work  for  possible  future 
projects,  for  a  computer  programme  which  estimates  the  accuracy  with  which 
the  parameters  of  a  hypothetical  satellite  orbit  could  be  determined  by  use  of 
ground  station  observations  of  one  or  more  specified  types.  Such  a  programme 
has  now  been  written  for  the  Pegasus  computer  and  is  described  in  this  Report. 

The  function  of  the  new  programme  is  most  easily  indioated  by  comparing 
the  programme  with  the  standard  R.A.E.  programme  for  orbit  determination  . 

The  two  main  differences  of  the  new  programme  are  the  following:  - 

(i)  it  is  not  confined  to  angular  observations,  but  can  deal  equally 
well  with  observations  of  range,  range  rate,  direction  cosines,  eto., 

(ii)  as  an  assessment  programme,  it  does  not  prooess  actual  observations 
of  an  aotual  satellite. 

The  significance  of  feature  (ii)  is  that  the  programme  is  much  simpler  to  use 
than  a  full  orbit  determination  programme  would  be,  and  it  requires  less 
computer  time. 

The  type  of  problem  to  which  the  new  programme  refers  may  be  summarised 
as  follows.  Nominal  values  of  K  unknown  parameters  of  a  hypothetical  satellite 
orbit  are  known,  as  are  also  the  locations  of  possible  tracking  stations.  The 
particular  quantities  which  each  station  can  measure  are  specified,  together 
with  the  nominal  accuracy  associated  with  each  quantity.  It  is  required  to 
estimate  the  accuracy  with  which  the  K  orbital  parameters  could  be  computed 
from  hypothetical  observations  over  a  given  number  of  days. 

The  difference  between  the  above  requirement  and  that  satisfied  by  a 
full  orbit  determination  programme  may  be  stated,  statistically,  rather  simply. 
All  observations  being  subject  to  random  errors,  assumed  independent  and 
normally  distributed,  we  may  consider  the  K-dimensional  variate  of  errors  in 
orbital  parameters  computed  from  a  given  set  of  observations.  An  orbit 
determination  programme  estimates,  by  maximum  likelihood,  the  first  and  second 
moments  of  this  variate.  An  assessment  programme,  on  the  other  hand,  estimates 
second  moments  only. 

The  pregramme  is  described  in  detail  in  Section  2  of  this  Report.  An 
application  to  an  assessment  problem  of  recent  interest,  involving  six-hour 
orbits,  is  described  in  Section  3. 

One  important  point  must  be  made.  Although  the  programme  does  not 
refer  to  aotual  observations,  the  computer  must  be  told  the  times  at  vhich  it 


k 


is  supposed  that  observations  could  be  made.  Problems,  such  as  (a)  finding 
the  periods  during  which  the  satellite  would  be  above  the  horizon  at  each 
station  and  (b)  deciding  how  many  independent  observations  to  assume in  such 
periods,  are  not  catered  for. 


2  THE  COMPUTER  PROGRAMME 
2,1  A  brief  description 


The  programme  in  question  is  entitled  "00,13.29  Error  Covariance  Estima¬ 
tion  for  Hypothetical  Orbit  Determination",  The  quantities  for  which 
covariance  estimates  are  made  are  certain  of  the  parameters  in  the  Merson 
model  used  at  R.A.E.  as  the  basis  for  satellite  orbit  determination  ,  As  in 
Ref,1  the  model  parameters  are,  using  standard  notation,  aQ,  e^,  i^,  0^,  tQ 

and  such  additional  parameters  e i.,  £!.,  u).  and  n.  as  are  required;  the  j 

0  J  J  J  0 

suffixed  parameters  allow  polynomial  contributions  to  orbital  elements  at  time  t, 

e»g,  2  e.(t  -  t  )J  to  eccentricity  e  and  2  n.("t;  -  t  )J  to  mean  motion  n 
JO  J  o 

(no  is  given  from  a^  by  Kepler* s  third  law).  Some  of  the  parameters  are 
regarded  as  known  quantities  and  are  such  that  in  a  genuine  orbit  determination 
they  would  be  held  fixed  at  pre-assigned  values.  The  remainder  are  the 
’unknown'  parameters  which,  starting  from  nominal  values,  it  would  be  the 
business  of  an  orbit  determination  programme  to  improve  by  use  of  observations. 
Programme  00,13,29  assesses  the  accuracy  with  which  the  unknown  parameters, 

K  in  number,  could  be  computed  in  such  determination.  The  distinction  between 
the  two  types  of  parameter  in  the  orbital  model  is  one  that  has  been  made  in 
the  working  version  of  the  R.A.E.  programme  for  orbital  determination  using 
angular  observations  only,  but  it  was  not  present  in  the  original  version  of 
that  programme  described  in  Ref.l. 


Input  for  the  programme  consists  of  the  assumed  values  of  all  parameters 
of  the  hypothetical  orbit  (including  the  fixed  ones)  and  of  information 
relating  to  a  number,  b,  of  possible  bursts  of  observational  data.  A  ’burst’ 
is  a  sequence  of  observations  from  the  same  station,  uniformly  separated  in 
time;  it  consists  of  (a)  station  position  data,  (b)  type  number,  indicating 
the  nature  of  the  quantities  which  could  be  measured  at  the  station,  (c) 
information  as  to  the  number  of  observations  in  the  burst  and  the  set  of  uni¬ 
formly  separated  times,  and  (d)  the  standard  deviations  of  errors  in  the 
quantities  measured.  As  remarked  in  Section  1  (and  as  explained  in  connection 
with  another  programme  on  p.1A  of  Ref .2)  no  values  of  hypothetical  observed 
quantities  themselves  have  to  be  given. 
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Output  consists  basically  of  the  standard  deviations  associated  with 

hypothetical  determination  of  all  the  unknown  parameters.  Frinting  is  in 

2 

the  same  standard  format  as  is  used  for  input  of  assumed  values,  zeros  being 

allocated  to  the  fixed  parameters.  Additional,  but  optional,  output  consists 

of  the  scaled  covariance  matrix  for  the  unknown  parameters.  This  is  a  K  x  K 

matrix  and  is  such  that  the  square  roots  of  the  diagonal  elements,  when 

de-scaled,  are  precisely  the  non-zero  standard  deviations  already  referred  to. 

The  scaling  is  irrelevant  since  the  only  likely  application  of  the  covariance 

2 

matrix  is  for  re-input  into 'the  Pegasus  programme  for  estimating  - 
errors  in  the  ephemeris  of  a  satellite  (see,  in  particular,  Section  2.5  of 
Ref. 2).:  the  scaling  is  then  automatically  allowed  for. 

2.2  General  theory 

Let  the  ’unknown*  parameters  of  the  satellite  orbit  be  E^.,  for 
k  =  1,2,...,K,  The  notation  E^.  will  also  be  used  for  nominal  values  or 
initial  estimates  of  these  parameters. 

Each  of  the  b  bursts  of  hypothetical  observations  is  associated  with  a 

particular  station.  Let  designate  the  station  associated  with  the  ith 

burst,  it  being  understood  that  a  station  receives  more  than  one  if  it  is 

responsible  for  more  than  one  burst.  Let  N.  denote  the  number  of  observations 

in  the  ith  burst  and  take  the  times  of  these  observations  to  be  t.  ,  t.  +  tt., 

t.  +  2t.,*..,  t.  +  (N.  -  1)  tt.  ,  where  t.  and  t.  are  known. 

1  111  1  1  1 

If  there  were  a  set  of  actual  observations  at  the  given  times,  they 
could  be  used  to  improve  the  orbital  parameters,  say  from  E,  to  E’ .  To  give 
the  theory  of  -the  assessment  problem  it  is  convenient  to  refer  to  the  hypo¬ 
thetical  observations  as  if  they  were  a  real  set.  This  set  will  include 
observational  errors  which  lead  directly  to  errors  in  the  improved  parameters 
E^..  However,  the  accuracy  assessed  for  the  E^.  and  given  by  standard 
deviations  and  covariance  matrix,  is  independent  of  the  actual  observations 
(and  errors)  and  depends  only  on  the  a  priori  accuracy  estimates  associated 
with  them.  It  is  assumed  here  that  the  E^  and  the  E*_  are  fairly  close,  so 
that  the  orbital  improvement  is  complete  after  one  iteration;  this  is 
equivalent  to  a  linearity  assumption  as  will  be  seen  in  due  course. 

Yle  suppose  then  that  we  have  a  set  of  observed  quantities.  Throughout 

the  ith  burst  let  the  s^  quantities  measured  be  denoted  by  0^,  where 

q  =  1,2,..,,  s^;  for  the  jth  observation  of  the  burst  we  suppose  that  0^  takes 

the  value  0  . this  occurs  at  time  t.  where  t. .  =  t.  +  (j  —  1)  m. .  With 
qij  *  ij’  10  1  1 

3^  =  2,  for  example ,  0^  might  be  right  ascension  and  0^  declination;  with 
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=  6,  0^  might  he  range,  02  and  direction  cosines  and  Q^,  0,.  and  0g  the 
time  rates  of  change  of  these  three  quantities.  The  possible  interpretations 
of  the  0  ore  listed  in  Section  2.3,  which  gives  formulae  connected  with 
observations  of  every  type  covered  by  the  programme. 

Although  the  three  suffices  q,  i  and  j  have  been  introduced  for 

precision,  it  will  be  convenient  in  most  of  the  following  analysis  to  omit 

them.  Y/e  concentrate  on  a  single  typical  observation  0  (=  ©  .  .)  made  by  a 

qm  j 

station  S»  We  suppose  that  0  contains  an  error  which  derives  from  a  normal 
distribution  having  zero  mean  and  the  particular  a  priori  standard  deviation 
associated  with  S. 

The  actual  observation  0  may  be  compared  with  a  theoretical  or  'computed* 
value,  0c,  which  is  a  function  of  orbital  parameters,  station  position  and 
time,  say 

eo  ■  s> 

If  limitations  of  the  orbital  dynamic  model  are  neglected,  ©c  is  a  known 

function  of  K  +  4  parameters,  since  E^,  S  and  t  account  for  K,  3  and  1 

respectively.  However,  we  make  two  further  assumptions:  that  the  co-ordinates 

of  the  given  station  are  known  without  error;  and  that,  in  contrast  with  the 

1 

basio  R.A.E.  orbit  determination  programme  as  originally  described  ,  time  also 
is  error-free.  Then  only  the  K  orbital  parameters  are  in  question  and  we 
now  write 

6o  *  60<V* 


The  comparison  of  0  with  ©c  provides  the  residual  R,  where 

R  *  0  -  0  (E.  ) . 

o'  k 

We  are  postulating  the  existence  of  a  complete  set  of  observations  from 
which  an  orbit  improvement,  of  to  E* ,  is  to  be  carried  out.  Writing 


it  follows  from  Taylor's  theorem  that,  to  first  order. 


R«  =  R  -  >  ■—  AE,  . 

/  r  0Ek  k 


Second  order  terms  may  be  nelgected  since  the  AE,  are  small;  this  is  the 

**  n. 

linearity  assumption  referred  to  earlier. 


Equation  (l)  can  be  set  up  for  each  0  ^  with  its  associated  <7^  (o'  is 
independent  of  suffix  j).  The  AE^.  are  then  found,  by  the  method  of  least 
squares,  from:- 


(R*  /c X  S 
qxj  qi 


=  minimum. 
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This  equation  is  solved  by  differentiating  with  respect  to  each  AE^.,  using 
equation  (1).  The  resulting 'equations  of  condition*  are  best  expressed  in 
matrix  form,  following  the  notation  of  Ref ,1  as  far  as  possible. 

Let  Y,  Z  and  M  be  the  two  column  vectors  and  matrix  defined  by 


Y 


(R  ,  ./cr  . ) , 
v  qij  qx  * 


(AE,)  and  M 


where  we  regard  qij  as  a  single  (row)  suffix.  The  matrix  equation  of 
condition  is  then  (T  denoting  transposition) 


MT(Y  -HZ)  =  0 


and  the  solution  for  Z  is  given  by 

z  =  (mt  m)”1  mt  y. 

To  obtain  the  required  covariance  matrix  of  the  computed  orbital 
parameters  -  and  this  is  just  cov  Z  -  we  introduce  the  notation  Y  for  the 
column  vector  which  would  be  derived  from  error-free  observations  and  Z  for 
the  corresponding  AE^.,  the  E^  in  this  case  being  the  true  value  of  the  para¬ 
meters.  Then  Y  and  Z  are  the  means  of  populations  of  all  possible  Y  and  Z 
based  on  the  distribution  of  errors  in  observations.  So,  if  tt  denotes 
expectation, 

cov  Z  =  At  ((Z  -  Z)(ZT  -  ZT)] 

=  (mt  m)”*1  mt  n\(Y  -  y)(yt  -  yt)|  m(mt  m)"'1 

T 

(since  M  M  is  a  symmetric  matrix) 

=  (MT  M)**1  MT  cov  Y  M(MT  M)“\ 
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But  the  elements  of  Y  are  so  weighted  that  the  standard  deviation  of  each  is  1, 
and  we  assume  that  these  elements  are  uncorrelated.  Thus  cov  Y  is  the  unit 
matrix.  Hence 

cov  Z  =  (MT  M)  ' .  (2) 


Formula  (2)  is  the  basis  of  the  Pegasus  computer  programme.  As  desired, 

it  is  independent  of  any  actual  observations  (it  does  not  contain  Y).  The 

only  quantities  required  are  partial  derivatives  and  a  priori  standard 

deviations  of  error.  The  difference  between  this  assessment  formula  and  the 

formula  for  cov  Z  which  would  be  used  for  the  analysis  of  actual  satellite 

2 

observations  is  that  the  latter  contains  an  additional  factor  (the  e  of  Ref.l). 
This  (scalar)  factor  would  be  derived  from  final  residuals  and  would  enable  a 
priori  error  estimates  for  the  observations  to  be  converted  into  a  posteriori 
estimate s. 


Formulae  for  observation 


The  various  interpretations  for  the  observed  quantities  0  must  now  be 

listed.  Formulae  involved  in  programming  the  derivatives  *rrr-  will  be  given. 

dEk 

Depending  on  the  type  of  observing  equipment  used,  observed  quantities 
fall  into  six  basic  categories:  range,  angles,  direction  cosines,  range  rate, 
angle  rates  and  direction  cosine  rates.  Although  angles  can  be  obtained  at 
once  from  direction  cosines,  these  two  categories  must  be  treated  separately 
because  of  the  weighting  question.  Thus  if  a  station  observes  direction 
cosines  a  fixed  standard  deviation  is  assumed  for  all  such  observations;  the 


standard  deviations  of  the  two  angles,  to  which  a  given  pair  of  direction 
cosines  could  be  converted,  would  depend  on  the  particular  observed  direction 
and,  in  addition,  the  angular  errors  would  in  general  be  correlated. 


Range  measurement  gives  a  single  quantity,  denoted  here  by  p.  Angle 
measurement  yields  two  quantities  for  specifying  a  line  of  sight.  These  may 
be  right  ascension  and  declination,  a  and  6,  or  azimuth  and. elevation. 

However,  no  distinction  is  necessary  as  will  now  be  demonstrated.  Suppose 
a  and  8  are  the  fundamental  (independent)  quantities;  then  5  is  measurable  to 
the  given  fixed  accuracy,  Cg,  associated  with  the  station,  but  a  only  to  the 
variable  accuracy,  c r  ,  given  by 


c  •“  crp  sec  8  . 
a  o 


This  is  equivalent  to  a  circular  normal  distribution  for  errors  in  any  plane 
perpendicular  to  the  line  of  sight  and  so  equivalent,  again,  to  the  following 


9 


for  azimuth  and  elevation:  c(el)  fixed,  cr(az)  =  cr(el)  x  sec(el).  Thus 
azimuth  and  elevation  are  effectively  fundamental  quantities  also. 

With  directiom  cosines  there  is  a  similar  situation,  but  only  so  long 
as  we  assume  that  these  cosines  are  always  measured  with  respect  to  axes  in 
the  horizon  plane.  This  is  normally  the  case  in  practice,  for  example  with 
Mini track  interferometer  observations.  Hence  we  may  conveniently  suppose  the 
fundamental  quantities  to  be  Z  and  m,  the  direction  cosines  of  a  line  of  sight 

with  respect  to  ground  plane  north  and  east  axes  respectively.  The  case  of 

axes  in  some  plane  other  than  the  horizon  (ground)  plane  is  not  covered  by  the 
programme. 

The  remaining  quantities  which  may  be  observed  fall  into  the  three  rate 
categories.  They  are  denoted  here  by  p,  a  and  8, and  Z  and  m. 

The  Hype*  of  an  observation  is  determined  by  the  categories  of  "the 
quantities  which  are  included  in  the  observation.  If  it  is  assumed  that 
angles  and  direction  cosines  are  never  included  in  the  same  observation,  there 
remain  27  possible  types.  Of  these,  15  have  been  specifically  catered  for 
by  the  programme.  The  type  no.  of  each  of  these  is  listed  in  the  table  below, 

together  with  the  quantities  covered  -  i.e.  interpretations  of  0  -  for  eaoh 

type.  It  may  occasionally  be  necessary  to  use  the  programme  with  observations 
of  a  type  not  covered  by  the  table.  Suppose,  for  example,  that  j5  Z  and  m 

are  measured  (by  a  combination  of  Doppler  and  interferometer).  Then  the 

• 

observation  may  be  regarded  as  of  type  15  with  p  Z  and  m  measured  to  very  low 
accuracy  (say  <Xp  =  ICr  metres  and  =  or.  =  1000  (sec)  ),  All  the  12  types 
not  specifically  covered  in  the  programme  may  be  dealt  with  in  this  way. 


Observation  types 


Type  no. 

Quantities  observed 

Type  no. 

Quantities  observed 

1 

P 

8 

p,-e,m 

2 

a, 8 

9 

•  .  ? 

P,a,8 

3 

Z,  m 

10 

p,Z,m 

4 

11 

P,P 

5 

a,  6 

12 

♦ 

a?5  ?a,6 

6 

Z,  m 

13 

•  • 

•d,m,d,m 

7 

p  ,a,6 

14 

p,a,8,p,a,S 

15 

*  *  • 

p,  Z,m, p,Zfm 

Turning  our  attention  to  the  formulae  for  partial  derivatives,  we  let 
r,  a,  8  denote  geocentric  spherical  polar  co-ordinates,  corresponding  to  the 
topocentric  p,  a,  6.  Then  if  X,  Y,  Z  be  station  co-ordinates  in  the 


normal  geocentric  system  of  axes, 

p  oos  8  cos  a 

p  cos  6  sin  a 

and  p  sin  5 


=  r  cos  8  cos  a  ~  X, 
=  r  cos  8  sin  a  -  Y 


=  r  sin  8 


The  partial  derivatives  of  r,  a  and  8  depend  on  the  Merson  model  for 
satellite  motion  and  are  given  in  full  elsewhere  .  Assuming  these  derivatives 
to  be  known,  we  require  formulae  for  the  derivatives  of  observations  6  in  the 
six  basic  categories.  We  start  with  the  first  two  categories,  formulae  for 
the  derivatives  of  p,  a  and  8  being  obtainable  from  the  equations  just  given. 
For  generality  we  take  differentials  of  both  sides,  employing  matrix  notation. 
It  is  interesting  to  note  that  time  derivatives  can  most  sinqily  be  dealt  with 
by  use  of  moving  axes,  fixed  in  the  earth,  which  instantaneously  coincide  with 
the  normal  inertial  axes  (x  axis  towards  the  vernal  equinox  and  z  axis  towards 
the  north  pole).  In  this  case  AX  =  AY  =  AZ  =  0  and  if  GOg  is  the  angular 
velocity  of  the  earth  we  got 


-  cos  a  0  -  sin  a\  /  -  cos  8  0  sin  8 


-  sin  a  0  cos  a  (  sin  8  0  cos  8  j  p  cos  8(Aa-00j,  At) 

0  1  0  /  \  0  1  0  /  \  P  A8 


-  oos  a  0 


sin  a 


a\  /  -  cos  6  0  sin  6\  / 


cos  a 


0  /  \ 


sin  8  0  cos  8  J  i  r  cos  S(Aa-0Og  At)  •  (3) 

/  / 

0  1  0  /  \  r  A  8  / 


When  the  standard  Ii.A.E.  orbital  determination  programmes  were  written 
time  derivatives  were  needed  for  correction  of  observation  times  and  equation 
(3)  was  developed  with  the  At  terms  kept  in.  For  the  present  programme, 
despite  the  fact  that  observations  of  time  differentiated  quantities  are 
allowed  for,  formulae  for  time  derivatives  have  not  been  used.  Dropping  the 
At  terms  and  inserting  the  matrices  on  the  left-hand  side  of  equation  (3)  we 
get 


|  p  cos  8  Aa 
\  PAS  , 


=  UJ 


1  “2  “3  / 


r  cos  8  Aa 


For  the  matrix  of  partial  derivatives  with  respect  to  the  we  now  have,  at 
once. 


It  is  convenient  to  leave  the  factor  cos  8  with  the  a  derivatives  since,  as 
already  mentioned,  the  standard  deviation  of  each  measured  a  includes  a 

*•1  3 

sec  8  factor,  so  that  elements  cv  cos  8  -r=—  are  required  in  the  matrix  M  of 

8  dE, 

Section  2.2. 

The  partial  derivatives  of  t  and  m,  observations  in  the  third  basic 
category,  can  be  obtained  from  those  of  a  and  8.  Let  (3  be  the  astronimical 
latitude  of  the  observing  station,  XQ  its  longitude  (east)  and  S  the  sidereal 
time.  Define  X  by 

X  =  X  +  S  . 
o 
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Then 


-  sin  p  cos  X  -  sin  {3  sin  X  cos  p\  j  cos  8  cos  , 

—  sin  X  cos  X  0  1 1  cos  6  sin  a. 

oos  p  cos  X  cos  P  sin  X  sin  p/  \  sin  6 


where  the  third  direction  coane,  n,  is  introduced  here  only  to  give  a 
complete  expression.  The  3x3  matrix  corresponds  to  a  change  to  left-handed 

axes  but  this  is  of  no  significance  and  in  any  case  we  are  only  concerned  with 
the  first  two  rows  of  the  matrix. 


On  taking  differentials, 


(“\ 

Am 

\  An/ 


p  cos  X 

-  sin  p  sin  X 

cos  (3\ 

/-  sin  a 

-  sin  8  cos 

sin  X 

cos  X 

0 

I  cos  a 

-  sin  8  sin 

p  cos  X 

cos  p  sin  X 

sin  p/ 

\  0 

cos  8 

+  <x )  /  sin  (3  sin  X  -  sin  (3  cos  X  0\.  j cos  6  cos  At, 


-  oos  X 


-  sin  X 


cos  6  sin  a 


\-cos  (3  sin  X  cos  (3  cos  X  o/\  sin  6 

Dropping  the  terms  involving  n  or  At  and  introducing  the  partial  derivative 
notation,  we  have 


d£ 

dE 

dm 

dE 


1 


bl 

dEK 


1 


dm 

dE. 


K 


-sin  p  cos  X  -sin  p  sin  X  cosp\^  y-sina 
i  -  sin  X  cos  X  0 


.  c  \  /  c  da 

•sxno  cos  a\  ycos  o 


-sin  8  sin  a 
cos  6 


d8 

dE 


1 


cos  8 


da 

dE 


K 


db 

dE 

(5) 


K 
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Finally  we  consider  derivatives  of  quantities  from  the  last  three 
categories,  those  which  involve  rates  of  change  with  time.  It  is  clear  that 
formulae  for  second  derivatives  are  involved  since,  for  example, 

iL  -  iie_  .  -ife. . 

dE1  dE^dt  “  dtdE1 

It  is  not  difficult  to  extend  the  formulae  of  Merson  to  cover  the  required 
second  derivatives  of  geocentric  co-ordinates.  After  doing  so,  however,  it 
is  still  necessary  to  differentiate  the  right-hand  sides  of  equations  (4)  and 
(5)  to  get  the  second  derivatives  of  observed  quantities.  It  was  therefore 
decided  that  second  derivatives  should  be  computed  numerically,  using  the 
available  sub-routines  for  first  derivatives. 

Thus  for  an  observation  0  we  use  the  approximation 

^At  ,  (6) 

where  At  is  some  suitable  fixed  time  increment.  Unless  told  otherwise  the 
programme  takes  At  =  1  sec;  for  an  orbit  of  long  periodic  time,  a  larger  value 
of  At  would  be  more  appropriate. 

No  difficulty  arises  with  weighting  factors  for  observed  quantities  in 
the  fourth  and  sixth  categories.  For  Ihe  fifth  category,  however,  there  is 
some  uncertainty  as  to  the  appropriate  factor  to  associate  with  a.  The 
natural  factor  is  again  sec  8  on  the  basis  that 

o’.  =  crJ  sec  8 

a  6 

where  o*£  is  a  given  fixed  standard  deviation.  But  this  relation  would  be 
derived  from  cos  8  <x(a)  =  <x(§)  and  it  could  very  well  be  argued  that  one 
should  preferably  start  from  cr(cos  8.  a)  =  c(6).  The  difficulty  is  associated 
with  the  fact  that  radars  do  not  normally  make  measurements  of  angle  rates 
independently  of  measurements  of  the  angles  themselves.  (We  contrast  the 
situation  with  direction  cosines  when  $,  and  m  may  be  regarded  as  essentially 
uncorrelated  with  -6  and  m.)  Thus  fifth  category  observations  must  be  treated 
warily.  For  completeness,  the  programme  permits  their  use  and  does  so ■  by 
computing,'  instead  of  (6), 


ae(Ek,t  +  At)  ae(Ek,t) 


(7) 
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00  s 


^cos  6(t  +  At) 


da(Ek,t  +  At) 

3\ 


-  cos  8(t) 


3a(Ek,t) 


2*4  Operating  instructions 

In  the  rest  of  Section  2  it  is  assumed  that  the  reader  is  familiar  with 
the  use  of  the  Pegasus  digital  computer.  We  first  give  the  simple  operating 
instructions. 

(i)  Depress  hand-switch  0  and  the  hoot  key;  depress  switches  1,  2 
and/or  3  as  required  (see  below);  clear  the  others. 

Place  the  binary  programme  tape,  00.13*29, in  the  main  tape  reader. 

START  and  RUN  -  the  tape  is  read  and  a  77  (z)  stop  reached, 

(ii)  Place  the  parameter  tape  (see  Section  2.5)  in  the  main  reader. 

RUN  (or  Start  and  Run)  -  the  tape  is  read  in,  computing  proceeds 
and,  in  due  course,  results  are  punched  (see  Section  2.6). 

Computer  action  at  the  completion  of  a  run  (i.e.  after  punching  results 
for  a  given  parameter  tape)  depends  on  whether  hand-switch  3  is  set.  If  it  is, 
further  parameters  are  immediately  read  in  without  need  for  intervention;  the 
machine  can  be  left  unattented  while  many  sets  of  parameters,  all  on  one  tape, 
are  dealt  with.  If  switch  3  is  clear  the  computer  comes  to  a  77  stop  in  0.4, 
permitting  a  change  of  parameter  tape  if  further  runs  are  to  be  made.  The 
next  parameters  may  be  dealt  with  by  repeating  instructions  (ii). 

The  interpretation  of  hand- switch  2  is  explained  in  Section  2.5,  and 
of  hand-switch  1  in  Section  2.6.  It  will  be  found  that  switch  2  is  normally 
kept  clear. 

The  computer  time  required  to  run  the  programme  may  be  forecast  roughly 
as  follows;  for  each  observation  allow  20  seconds  if  its  type  does  not  include 
any  rate  measurement;  if  rate  measurement  is  included  allow  40  seconds. 

Two  amendments  to  the  programme  may  be  made,  if  desired,  and  fed  into 
the  computer  after  (i)  above.  They  are  referred  to  in  Sections  2.6  and  3.6 
respectively. 

The  value  of  At  (see  Section  2.3)  may  be  changed  from  1  second  to  any 
integral  multiple,  p,  of  10  |_isec.  Two  actions  are  necessary:-  (a)  the 
programme  must  be  amended  by  feeding  in  the  tape 

T  266.2 

+  P 

Z 
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(b)  the  standard  deviations  of  all  rate  measurements  must  be  multiplied  by 
p/100000  before  they  are  punched  on  the  parameter  tape  (see  Section  2.5). 

2.5  Parameter  tape 

The  print  out  of  an  illustrative  parameter  tape  is  given  below. 

The  numbers  on  "the  right  refer  to  blocks  of  punching,  separated  by  black  tape, 
and  are  introduced  here  for  reference  in  the  commentary  which  follows. 


T233-7 

(i) 

+2 

(no.  of  bursts) 

J300*0 

(ii) 

16730 

(semi-major  axis) 

0*0003 

(eccentricity) 

80 

(inclination) 

255 

(RA  of  node) 

0 

(argument  of  perig 
(time  at  node) 

1964  Jan  1  000 

0 

500 

* 

T227-0 

(iii) 

+15 

-454159256 

(type) 

00 

+427230370 

W 
,  (z)  . 

-133904980 

-0*211324797 

(sin  pj 

+0*977415893 

(cos  3) 

+0*685182992 

(sin 

-0-728370968 

(cos  X  ) 

L 

0 

1964  Jan  1  14  13  0 

11*5 

11 

(iv) 

+50 

+0*000175 

+1 

+0*0000077 

(v) 

1964  Jan  1 

11*5 

11 

+50 

+0*000175 

+1 

+0*0000077 


22  45  0 


(vi) 

(vii) 

(viii) 
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(i)  The  number,  b,  of  bursts  of  data  must  first  be  set  in  storage 
location  233.7;  here  b  =  2. 

(ii)  The  programme  is  entered  at  address  300.0  and  a  block  of  orbital 
parameters  is  at  once  read  in.  Quantities  on  consecutive  lines  are  aQ,  eQ> 
iQ,  o)Q  and  tQ,  the  0  and  500  being  (arbitrary)  dummy  numbers.  Polynomial 
coefficients  may  be  incorporated  in  the  parameter  format,  e.g.  e^  after  e  , 

n,j  and  n^  after  the  dummy  500.  For  complete  details  of  this  standard  parameter 
format  and  punching  rules.  Section  2.1+  of  Ref. 2  should  be  consulted. 

Ref. 2  makes  a  distinction  between  parameters  which  are  'btrict  E  type" 
and  those  which  are  not.  This  distinction  is  the  same  as  the  one  made,  in 
Section  2.1  of  the  present  paper,  between  parameters  which  are  'unknown'  and 
'fixed'  respectively.  In  the  normal  use  of  programme  00.13.29  there  will  be 

no  fixeS  parameters.  The  programme  can  still  be  used  when  there  are  one  or 

more,  however,  as  follows:- 

(a)  the  parameter  type  must  contain  an  additional  section;  between 
Sections  (ii)  and  (iii)  (after  "the  block  of  orbital  parameters)  must 
appear 


+0.uvwxyz 

L 

(b)  hand-switch  2  must  be  set  (enabling  the  new  section  of  tape  to  be 
read) . 

The  interpretation  of  u,...,z  is  as  in  Ref. 2,  but  it  is  noted  that  the 
additional  punching  occurs  at  different  parts  of  the  parameter  tape  in  the 
two  programmes. 


(iii)  Station  data  is  always  stored  in  block  227  of  the  computer. 

The  first  quantity  is  type  number,  15  -  designating  p,  Z,  m,  |3,  Z  and  m  - 
in  the  case  of  the  illustrative  example.  The  next  three  quantities  are 
station  co-ordinates  X,  Y  and  Z  in  cm.  The  remaining  four  are  sin  (3,  cos  (3, 
sin  XQ  and  cos  where  (3  and  Xq  are  latitude  and  longitude  (east). 

(iv)  The  times  of  the  observations  in  the  first  burst  of  data  are 
specified  by  punching  t^  (see  Section  2.2),  (in  minutes)  and  itself. 


(v)  Standard  deviations  for  errors  in  the  quantities  observed  are 


specified,  viz.  Cp,  cr ^ 


(=  cr  ),  cr.  and  cr*  (*>  cr. )  for  the  example  given.  Units 
v  m  *  p  Z  nr 

for  the  six  categories  of  observed  quantities  are:  for  p,  metres;  a  and  6, 

(5, metres/sec;  a  and  6,  seconds  of  arc/seo;  and 


seconds  of  arc;  Z  and  m. 
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•6  and  m,  -/sec;  the  standard  deviations  must  be  punched  in  the  order  of  the 
categories  -which  are  included  and  they  must  be  followed  by  an  asterisk.  (They 
are  read  in  by  the  Pegasus  Matrix  Interpretive  Scheme.) 

(vi)  The  computer  looks  for  station  data  for  the  second  burst,  similar 
to  data  under  (iii).  If  there  is  no  change,  either  in  observation  type  or 
co-ordinates,  it  is  only  necessary  to  punch  the  L.  This  is  a  return  link 
from  initial  orders,  used  as  a  sub-routine,  to  the  programme. 

(vii)  Information  for  the  second  burst,  -  tg^i  ^  and  Ng  -  is  punched 
as  under  (iv) . 

(viii)  The  standard  deviations  for  the  second  burst,  even  if  identical 
with  those  for  the  first,  must  be  punched  in  full. 

NB  Positive  signs  have  been  omitted  in  the  punching  of  Sections  (ii),  (iv) 
and  (vii),  where  they  are  optional.  In  Sections  (i),  (iii),  (v)  and  (viii) 
they  are  compulsory. 

2.6  Output 

The  regular  output  consists  of  a  block  of  parameter  standard  deviations, 
each  denoting  the  accuracy  with  which  the  corresponding  orbital  parameter  could 
be  computed  from  the  given  observational  data.  Thus  for  the  illustrative 
parameters  and  data  used  in  the  previous  section,  the  output  was  as  follows: 


+0-0478 
+0-000003 
+0-0008 
+0-0007 
+0-8461 
0  00  00-334 
+0- 0000 
+0-00 


(cra  in  km) 

(%) 

(o"i  in  degrees) 

(oq  in  degrees) 

(cr^  in  degrees) 

(o"^o  in  h  m  s) 

(associated  with  the 
dummy  0  and  500  of  input). 


If  any  parameter  is  of  the  fixed  variety  the  appropriate  standard  deviation 
is  printed  as  zero.  If  required,  each  quantity  can  be  printed  to  three 
further  decimal  places  by  feeding  in  an  amendment  to  the  basic  programme. 

If  the  computer  is  operated  with  hand-switch  1  clear  there  is  no  other 
printing  than  that  above.  If  this  switch  is  set,  however,  the  standard 
deviations  are  preceded  by  the  punching  of  the  complete  covariance  matrix, 
cov  Z,  from  which  they  are  derived.  There  is  little  point  in  giving  the 
complete  matrix  for  the  example  being  considered,  since  the  units  of  distance 
(a),  angle  and  time  are  not  -the  obvious  ones.  The  first  two  columns  are 
given  for  interest,  however,  in  the  standard  floating  point  (argument  exponent) 
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form*  The  element  1.073523  x  10  is,  of  course,  the  non-dimensional 

variance  cr  .  The  values  of  the  covariance  u  from  the  first  and  second 

T  P  P 

columns  differ  slightly  due  to  rounding  error.  From  cr  ,  u  and  cr  may  be 

a.  fits  © 

obtained  the  correlations  between  errors  in  semi-major  axis  and  eccentricity; 
its  magnitude  is  +0.23. 


1st  col.  +3*030766  -16 
+1*304584  -14 
-5*627968  -14 
+1-050549  -14 
-9*670146  -12 
-2-111723  -15 


2nd  col.  +1-304585  -14 
+1*073523  -11 
+1*771246  -13 
+1-227553  -11 
-5*823236  -10 
-9-643007  -14 


The  standard  deviations  are  valuable  in  that  they  provide  an  immediate 

insight  into  accuracy  attainable  with  a  given  set  of  observational  data. 

The  complete  covariance  matrix,  particularly  in  its  scaled  floating  point  form* 

adds  little  to  this  quick  picture.  Thus  even  when  the  matrix  is  included  in 

the  computer  output,  it  will  not  normally  be  sensible  to  print  the  tape.  The 

necessity  for  the  complete  matrix  arises  if  the  accuracy  (expressed  as  a 

standard  deviation)  with  which  some  function  of  the  computed  orbital  parameters 

is  known  is  to  be  assessed.  The  obvious  application  is  for  estimating  errors 

in  an  ephemeris  based  on  the  computed  parameters;  the  computer  programme 

2 

already  written  for  this  purpose  requires  the  cov  Z  output  tape  as  part  of 
its  input . 

3  AN  APPLICATION  OF  THE  PROGRAMME 
3.1  Background 

The  European  Launcher  Development  Organisation  (ELDO)  have  recently  been 
interested  in  the  accuracy  with  which  parameters  for  certain  orbits  could  be 
computed  from  tracking  data  generated  by  a  single  radar  set.  The  orbits 
have  been  of  the  circular  near-polar  type,  and  such  that  the  period,  as  seen 
from  the  earth,  is  a  simple  fraction  of  a  day.  Two  particular  cases  are  the 
*6  hour*  and  ’8  hour’  orbits,  when  the  fraction  is  ^  and  -g  respectively,  A 
’day*,  here,  must  be  taken  to  be  that  period,  almost  equal  to  a  sidereal  day 
but  with  allowance  for  satellite  orbit  precession,  after  which  the  track  of  the 
satellite  on  the  earth  repeats  itself. 


It  was  thought  to  be  sensible  that  the  first  application  of  the  new 
computer  programme  should  be  to  one  of  these  orbits  of  interest  to  ELDO.  The 
6  hour  orbit  was  chosen.  The  radar  was  taken  to  be  sited  in  the  north-east 
of  Australia. 


£ 
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For  the  particular  orbit  parameters  considered  (see  Section  3.2)  the 
satellite  could  be  tracked  twice  a  day,  first  when  going  south  and  second  when 
going  north  an  orbit  and  a  half  later.  It  was  decided  to  look  at  the  accuracy 
question  for  two  cases;  first,  from  the  point  of  view  of  quick  answers,  how 
accurate  an  orbit  could  be  obtained  from  data  over  the  first  (south-going) 
pass  or,  indeed,  over  less  than  this  full  pass;  seoond,  from  the  point  of  view 
of  knowledge  of  the  orbit  after  two  or  three  days,  what  would  be  the  accuracy 
when  data  from  several  passes  were  combined. 

3«2  Orbital  parameters,  radar  site  and  observation  type 

The  orbital  parameters  are  those  listed  in  Section  2.5  For  an  illustra¬ 
tive  parameter  tape.  For  a  6  hour  orbit  at  80°  inclination  the  effect  of 
precession  is  to  shorten  the  sidereal  day  by  0.24  min.  The  actual  orbital 
period  is  thus  5^58m57S#  to  which  corresponds  the  assumed  semi-major  axis  of 
16730  km.  Though  the  orbit  is  nominally  circular,  eccentricity  was  taken  to 
be  0,0003  to  avoid  any  possible  difficulty  with  e  factors.  The  argument 
of  perigee  is  of  course  arbitrary  and  was  taken  zero;  the  standard  deviation 
obtained  for  w  was  inevitably  going  to  be  large,  due  to  the  small  e,  but  this 
was  of  no  importance.  The  parameter  tQ,  time  at  the  node,  was  also  arbitrary 
and  taken  at  January  1.0  of  1964.  The  value  of  0  was,  of  course,  far  from 
arbitrary;  255°  was  derived  by  making  assumptions  about  the  geography  of 
launch  and  injection  and  by  adding  in  the  sidereal  time  at  tQ. 

For  a  high  orbit  there  is  no  need  to  admit  an  unknown  parameter  tQ 
represent  the  effect  of  air  drag;  • so  no  parameter,  other  than  the  basio  six, 
wa3  introduced.  For  a  similar  exercise  with  a  low  orbit  it  would  be 
necessary  to  introduce  a  parameter  n^, possibly  also  n the  value  of  would 
not  have  to  be  known  and  could  in  fact  be  set  to  zero,  the  important . thing 
being  the  presence  of  the  additional  parameter  in  the  model. 

The  position  of  the  radar,  in  Australia,  was  taken  to  be  at  about 
latitude  12°  S  and  longitude  137°E.  The  station  parameters  used  for  the 
programme  are  those  given  by  the  illustrative  parameter  tape  of  Section  2,5* 

The  radar  was  assumed  to  make  observations  of  type  5>  i.e.  to  measure 

the  following  six  quantities;  slant  range,  p;  direction  cosines  relative  to 

*  • 

ground  plane  axes,  &  and  m;  p;  £  and  m*  Accuracies  were  taken  to  be  given 

by  the  following  expressions  for  standard  deviations;  cr^  =  50m> 

cr.  -  cr  ~  1.75  x  ICf^;  or.  =  1  m/sec,  cj  =  cr.  =  7.7  x  10“  6/ sec. 

*0  m  P  Kj  in  ' 
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3.3  Observation  bursts 

For  •the  given  orbit,  as  looked  at  from  the  given  radar  site,  the 
satellite  would  appear  above  the  horizon  four  times  a  day,  twice  going  south 
and  twice  going  north.  One  of  the  north  passes  and  one  of  the  south  may  be 
ruled  out  at  once,  however,  on  the  ground  that  the  satellite  would  never  climb 
to  an  elevation  greater  than  10°.  Of  the  other  two  passes  the  satellite  would 
reach  about  50°  on  the  south-going  and  about  60°  on  the  north-going.  Taking 
10°  as  (an  arbitrary)  radar  horizon,  the  approximate  times  of  visibility  are 
as  follows: 


south  going 

14h13m 

to 

16h8”, 

north  going 

22V 

to 

*v. 

an  interval  of  1 1 5  minutes  in  eaoh 

case. 

The  se 

are  the  times  for  the  first 

day  (1964  January  1);  they  become  4m12S  earlier,  each  successive  day. 

It  was  decided,  as  stated  in  Section  3.1*  to  carry  out  the  computer  runs 
under  two  heads.  The  first  covered  the  use  of  data  from  just  the  first 
possible  pass,  starting  with  an  observation  at  14  13  ;  successive  observations 
were  at  intervals  of  5  minutes  (chosen  arbitrarily),  totals  up  to  a  maximum 
possible  of  24  observations  being  considered.  The  second  head  covered  runs 
using  more  than  one  burst  of  data;  each  burst  corresponded  to  one  pass  and  was 
divided  into  10  intervals  of  11.5  minutes,  giving  11  observations  over  the 
burst.  Under  the  second  head  runs  were  carried  out  for  one  burst  only,  two, 
three,  four  and  six;  the  final  case  corresponded  to  use  of  data  from  all 
passes  over  a  three  day  period. 

3.4  Results  for  data  using  first  pass  only 

Let  N  be  the  number  of  observations  used,  a  suffix  being  redundant  in 
the  case  of  a  single  pass.  The  observations  are  assumed  to  be  made  at 
5  minute  intervals  so  that  the  last  corresponds  to  a  time  5(N-l)  minutes  from 
the  (effective)  beginning  of  the  pass.  The  maximum  value  of  N  is  24.  For 
values  of  N  less  than  the  maximum  the  programme  shows  how  the  accuracy  of  orbit 
determination,  using  data  up  to  a  given  point  reached  in  the  pass,  varies  as 
that  point  advances  in  time. 

Results  were  obtained  for  a  number  of  values  of  N  from  5  to  24.  The 
corresponding  standard  deviation  estimates  cr^,  cr^,  o\ ,  c ecr^  and  cr-fc0  are 
plotted  against  N  in  Figs. 1-6.  The  reason  for  plotting  e  cr  ,  with  e  =  0.0003, 
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has  "been  implied  in  Section  3*2.  It  can  be  seen  clearly  if  it  is  observed 
that  standard  deviations  in  the  non- singular  orbital  elements  c(=e  cos  go)  and 
s(=esin  go)  are  given,  when  go  =  0,  by 

cr  =  o'  and  o  =  e  o'  . 
c  e  s  go 

It  is  of  interest  to  know  how  far  the  decrease  in  each  o  for  increasing 
N,  indicated  in  Figs. 1-6,  is  caused  by  the  fact  that  an  increasingly  long  arc 
of  the  orbit  becomes  available  and  how  far  by  the  mere  fact  that  there  are  more 
observations.  Since  the  standard  deviation  of  any  measured  quantity  is  in 
general  inversely  proportional  to  the  square  root  of  the  number  of  observations 
made,  we  can  compute  an  ’accuracy  per  unit  observation*  by  multiplying  each  c r 
by  the  appropriate  /H.  For  semi-major  axis.  Fig. 7  has  been  obtained  in  this 
way  from  Fig.1.  It  confirms  a  result  that  might  be  expected  intuitively, 
namely,  that  increasing  coverage  is  important  until  about  half  the  pass  has 
been  seen,  while  data  from  the  second  half  gives  small  return. 

Results  from  just  one  pass  -  less  than  two  hours  worth  of  data  -  are 
thus  surprisingly  good.  From  merely  24  observations  semi  major  axis  is  known 
to  about  250  metres.  Hence  orbital  period  is  known  to  about  0.5  sec.  This 
in  turn  means  that  the  suspected  drift  of  the  satellite  after  a  full  (modified 
sidereal)  day  is  known  to  about  2’  arc  or,  after  a  year,  to  about  12°.  It  is 
interesting  to  know  to  which  of  the  four  sorts  of  quantities  being  measured  - 
p,  Z  and  m,  p,  Z  and  m  -  the  good  accuracy  is  due.  To  throw  light  on  this 
point  the  following  repeats  of  the  24  point  run  were  carried  out,  with 
different  assumptions  as  to  the  observation  type:-  (a)  type  1  (p  only), 

(b)  type  3  (Z  and  m),  (c)  type  4  ((5  only),  (d)  type  6  (Z  and  m),  (e)  type  11 
(p  and  p),  (f)  type  13  (Z,  m,  Z  and  m) .  The  results  are  listed  in  Table  1. 

An  immediate  conclusion  is  that  measurement  of  Z  and  m  is  of  little  value  and 
adds  nothing  to  measurement  of  Z  and  m.  In  fact  to  "pull  their  weight"  the 
accuracy  of  angular  rate  observations  would  have  to  be  improved  by  a  factor  of 
about  20.  Further  features  of  Table  1  are  discussed  in  Section  3*6* 

3.5  Results  for  data  using  several  passes 

Results  were  obtained  in  a  similar  way  to  that  of  the  previous  section* 
Complete  data  over  each  of  several  passes  being  assumed,  the  number  of  passes 
used  was  increased  from  one  to  six.  To  restrict  the  use  of  computer  time  the 
number  of  observations  per  burst  was  cut  to  11  leaving,  even  then,  66  observa¬ 
tions  to  be  dealt  with  in  the  case  of  six  passes  analysed.  For  the  analysis 
of  the  first  pass  only,  of  course,  nothing  was  added  to  previous  knowledge. 


It  is  meaningless  to  give  the  results  in  graphical  form.  Instead,  they 
are  set  out  in  Table  2;  the  results  for  the  case  of  two  passes  have  been  given 
also  as  the  illustrative  output  in  Section  2.6.  The  table  agrees  well  with 
intuitive  notions.  Thus  as  soon  as  data  from  a  second  pass  are  included,  this 
pass  being  some  eight  hours  later  than  the  first  and  covering  a  very  different 
part  of  the  orbit,  the  accuracy  of  all  parameters  is  improved,  a  factor  of  about 
5  or  6  being  involved  in  each  case.  Addition  of  a  third  pass,  however,  affects 
the  accuracies  of  the  elements  in  quite  different  ways.  Eccentricity  and  the 
angular  elements  -  i,  fi  and  go  -  are  scarcely  improved  since  the  third  pass  merely 
repeats  the  first.  But  semi-major  axis  is  improved  by  a  factor  of  nearly  20, 
because  for  the  first  time  a  really  accurate  measure  of  orbital  period  becomes 
available.  The  effect  on  the  sixth  element,  tQ,  is  nearly  as  good.  The 
improvement  in  semi-major  axis  as  up  to  six  passes  are  included  fits  in  with 
the  view  that  the  standard  deviation  of  the  error  with  which  orbital  period 
may  be  measured  is  inversely  proportional  to  the  total  time  over  which 
observations  are  used.  A  standard  deviation  for  semi-major  axis  of  less 
than  1  metre  (after  3  days)  is,  of  course,  remarkably  good.  It  corresponds 
to  a  drift  in  satellite  position  of  no  more  than  0°.05  at  the  end  of  a  year. 

3.6  Further  discussion 

It  is  wished  to  draw  attention  to  three  points  which  arose  in  connection 
with  the  present  application  of  the  computer  programme  but  which  have  general 
significance. 

The  first  point  concerns  the  importance  of  the  non-diagonal  terms  of  the 
covariance  matrix,  cov  Z.  It  is  well  brought  out  by  reference  to  Table  1, 
introduced  in  Section  3.4.  It  may  be  wondered  why  it  is  that  use  of  obsepva- 
tions  from  all  of  the  four  relevant  categories  (p,  £  and  m,  p,  £  and  m)  yield 
so  much  better  accuracy  km)  for  semi-major  axis,  in  particular,  than 
observations  from  any  single  category  (15  km  at  best). 

If  C.j,  C2,  C^,  denote  the  covariance  matrix  for  observations  restricted 
to  each  of  the  four  relevant  categories  in  turn,  then  the  overall  covariance 
matrix  is  given  by 

(cov  Z)"1  =  C"1  +  C"1  +  C'1  +  C”1. 

Hence  if  non  diagonal  terms  of  the  matrices  were  neglected  we  should  expect, 
from  Table  1 , 

o-"2  cs  38"2  +  15~2  +  32~2  +  496 "2  =  (13)“2, 

& 

so  that  the  actual  cr  is  about  $0  times  better  than  expectation* 

a 
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The  explanation  lies  in  the  relevant  magnitude  of  the  non-diagonal  terms, 

i»e.  in  the  correlations  between  semi-major  axis  and  the  other  five  elements. 

These  correlations  are  large  for  all  the  presently  considered  matrices  C^,  C^, 

Cy  and  cov  Z.  For  cov  Z  (for  which  <ra  =  \  km)  the  correlation  between 

a  and  e  is  -0.75  and  between  a  and  t  -0.964. 

o 

The  second  point  also  relates  to  the  presence  of  large  correlation 
elements  in  the  covariance  matrices.  As  seen  from  equation  (2),  the  final 
operation  in  the  determination  of  cov  Z  is  the  inversion  of  the  matrix  (M  M). 

If  this  matrix  is  ill-conditioned  there  may  be  considerable  loss  of  accuraoy 
in  the  performance  of  the  operation. 

The  matrices  leading  to  the  results  listed  in  Table  1  were  all  rather 
badly  conditioned.  The  outstanding  case,  surprisingly,  was  for  the  first 

line  of  the  table,  corresponding  to  orbital  determination  from  range  data  only. 

T 

The  standard  deviations  given  for  this  case  had  to  be  found  by  inverting  (M  M) 
with  the  aid  of  the  Double  length  Matrix  Interpretive  Scheme  for  Pegasus. 

Use  of  the  normal  scheme  (incorporated  in  the  programme)  had  previously  led  to 
results  which  underestimated  all  the  standard  deviations  by  the  same  factor, 

0.42;  c r  ,  for  example,  had  been  estimated  at  1'6  km. 

One  reason  for  ill-conditioning,  in  the  case  under  consideration,  may  be 
seen  at  once.  The  sixth  element  in  the  orbital  model,  t  ,  is  a  time  which  is 
14  hours  before  the  beginning  of  the  2  hour  burst  of  data  used  to  determine 
the  elements.  The  situation  is  analogous  to  the  difficulty  encountered  if  one 
tries  to  fit  a  straight  line  y  =  mx  +  c  to  a  set  of  points  (x,y)  which  are  all 
bunched  well  to  one  side  of  the  y-axis. 

If  ill-conditioning  is  suspected  when  the  programme  is  to  be  used,  the 
remedy  is  to  feed  in  an  amendment  to  the  programme  such  that  the  matrix  (M  M) 

is  output  instead  of  (or  in  addition  to)  standard  deviations  and  cov  Z. 

This  matrix  may  then  be  inverted  by  means  of  a  short,  specially  written  programme 
which  uses  double-length  arithmetic  with  a  working  length  of  22  significant 
figures  instead  of  only  9. 

The  final  point  concerns  the  choice  of  dynamic  model.  Let  us  suppose 
that,  for  a  particular  orbit,  no  parameter  other  than  semi-major  axis,  say, 
is  of  interest.  It  is  tempting  then  to  argue:  "since  we  have  nominal 

(e.g.  launch)  values  of  e,  i,  fi,  00  and  tQ,  and  since  we  are  not  interested  in 

improving  these,  let  us  keep  them  fixed  and  see  how  accurately  we  can  determine 
the  one  parameter  we  want  to  know". 
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In  the  case  of  the  example  "being  considered  this  would  involve  use  of  a 

parameter  tape  containing  +0.100000  (see  Section  2.5);  values  for  cr  far 

a 

better  than  those  quoted  in  this  paper  would  be  obtained. 

The  snag,  of  course,  is  that  if  e  etc,  are  held  fixed  they  are  assumed 
known  exactly,  whereas  they  actually  will  contain  errors  which  react  on  a. 

This  reaction  will  be  important  unless  e  etc,  are  known  to  greater  accuracy 
than  that  to  which  they  could  be  computed  if  treated  as  unknowns.  Such 
knowledge  will  not  normally  exist,  though  it  has  been  assumed  as  a  justifica¬ 
tion  for  omitting  the  parameter  n^  (see  Section  3,’l)  in  the  present  exercise. 
Summarising  for  the  general  case,  of  which  the  present  example  is  typical: 
an  assessment  study  must  assume  the  necessity  to  estimate  all  parameters,  even 
though  the  answers  for  only  one  or  two  may  be  actually  looked  at. 

4  CONCLUSIONS 

A  programme  has  been  written  which  fills  an  important  gap  in  the  Pegasus 
computer  library  of  satellite  orbit  programmes.  It  should  be  useful  in 
assessment  work  on  the  accuracy  with  which  hypothetical  orbits  can  be  deter¬ 
mined  from  hypothetical  observations. 

A  first  application  of  the  programme  has  demonstrated  that  accurate 
determination  of  circular  near-polar  six-hour  orbits  can  be  made  using  radar 
data  from  a  single  station  in  Australia.  This  example  has,  at  the  same  time, 
focussed  attention  on  the  danger  which  is  inherent  in  the  operation  of 
inverting  ill-conditioned  matrices. 


Table  1 

ACCURACY  OF  ORBIT  DETERMINATION  FROM  SINGLE-PASS  DATA 


FROM  OBSERVATIONS  OF  VARIOUS  TYPES 


Observation  ^ 
type  a 


1  (p  only)  38 

3  (£  and  m)  15 

4  (p  only)  32 

6  (G  and  m)  496 

11  (p  and  |5)  21 

13  (G,m,G,m)  15 

5 


o\  (deg) 


0.0034  0.798 
0.0007  0.005 
0.0022  1.398 

0.0234  0.442 
0.0018  0.439 
0.0007  0.005 
4 


Oq  (deg) 

eorw  (deg) 

°t0  (sec) 

0.656 

0.121 

249 

0.004 

0.014 

73 

0.391 

0.075 

162 

0.616 

0.586 

2502 

0.352 

0.064 

137 

0.004 

0.014 

73 

0.002 

0.001 

1 

1 - 

Table  2 

ACCURACY  OF  ORBIT  DETERMINATION  FROM  MULTI-PASS  DATA 


passes 

used 


cr  (metres) 
a 

107cr 

e 

376 

164 

48 

33 

2,5 

30 

1.8 

25 

0.9 

20 

c (sec  ) 
o 


1.7 

0.3 

0.024 

0.021 

0.016 


f. 

V:  r  -a  - 
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FIG.  I  ACCURACY  OF  SEMI-MAJOR  AXIS  DETERMINATION  V. 
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FIG. 2  ACCURACY  OF  ECCENTRICITY  DETERMINATION  V. 
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FIG. 3  ACCURACY  OF  INCLINATION  DETERMINATION  V. 
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FIG.4  ACCURACY  OF  NODAL  RIGHT  ASCENSION  DETERMINATION  V. 
NUMBER  OF  OBSERVATIONS  5m  APART 
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FIG. 5  ACCURACY  OF  ARGUMENT  OF  PERIGEE  DETERMINATION  V. 
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FIG. 6  ACCURACY  OF  NODAL  EPOCH  DETERMINATION  V. 
NUMBER  OF  OBSERVATIONS  5m  APART 


8000 

tf'ajN" 

(metres) 

4000 


o 


DETACHABLE  ABSTRACT  CARDS 


from  observational  data  of  specified  type  and  assumed  accuracy#  An  from  observational  data  of  specified  type  and  assumed  accuracy.  An 
application  of  the  progranme  is  made  to  an  orbit  of  six  hours  period  application  of  the  programme  is  made  Xo  an  orbit  of  six  hours  period 
determined  from  radar  observations  at  a  single  station*  determined  from  radar  observations  at  a  single  station. 


